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Abstract 

In this note, we show how certain properties of Goldbeter's 1995 model 
for circadian oscillations can be proved mathematically, using techniques 
from the recently developed theory of monotone systems with inputs and 
outputs. The theory establishes global asymptotic stability, and in partic- 
ular no oscillations, if the rate of transcription is somewhat smaller than 
that assumed by Goldbeter. This stability persists even under arbitrary 
delays in the feedback loop. 



1 Introduction 

The molecular biology underlying the circadian rhythm in Drosophila is the 
focus of a large amount of both experimental and theoretical work. Goldbeter 
proposed a simple model for circadian oscillations in 0] (see also his book [5]). 
Although by now several more realistic models are available, in particular incor- 
porating other genes, this simpler model exhbits many realistic features, such as 
a 24-hour period. The key to the model is the inhibition of per gene transcrip- 
tion by its protein product PER, forming an autoregulatory negative feedback 
loop. 

In this note, we show how certain properties of the model can be proved 
mathematically, using techniques from the recently developed theory of mono- 
tone systems with inputs and outputs. The theory establishes global asymptotic 
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stability, and in particular no oscillations, if the rate of transcription is some- 
what smaller than that assumed by Goldbeter. This stability persists even under 
arbitrary delays in the negative feedback loop. On the other hand, a larger - 
but still smaller than Goldbeter's- strength, in the presence of delays, results 
in oscillations. 

The terminology and notations are as given in [21 EI ? and are not repeated 
here. 



2 The Model 

The model is as shown in Figure ^ PER protein is synthesized at a rate pro- 
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Figure 1: Goldbeter's Model 

portional to its mRNA concentration. Two phosphorylation sites are available, 
and constitutive phosphorylation and dephosphorylation occur with saturation 
dynamics, at maximum rate tij's and with Michaelis constants Ki. Doubly 
phosphorylated PER is degraded, also satisfying saturation dynamics (with pa- 
rameters Vd,kd), and it is translocated to the nucleus with rate constant k\. 
Nuclear PER inhibits transcription of the per gene, with a Hill-type reaction 
of cooperativity degree n and threshold constant Ki, and mRNA is produced, 
and translocated to the cytoplasm, at a rate determined by a constant v s . Ad- 
ditionally, there is saturated degradation of mRNA (constants v m and k m ). 
The equations for concentrations are as follows: 

M = v s K?/(K? + P%)-v m M/(k m +M) 

P = k s M-V 1 P /(K 1 +P Q ) + V 2 P 1 /(K 2 +P 1 ) 

Pi = V 1 P /(K 1 + P Q )-V 2 P 1 /(K 2 + P 1 )-V 3 P 1 /(K 3 + P 1 ) + V i P 2 /(K 4 + P 2 ) 

P 2 = V 3 P 1 /(K 3 +P 1 ) - V 4 P 2 /(K A + P 2 ) - k x P 2 + k 2 P N - v d P 2 /(k d +P 2 ) 

Pn - hP 2 -k 2 P N 



where the subscript i = 0, 1,2 in the concentration Pi indicates the degree of 
phosphorylation of PER protein, Pn is used to indicate the concentration of 
PER in the nucleus, and M indicates the concentration of per mRNA. The 
parameters (in suitable units \iM or are as in the following table: 
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Parameter 


Value 


Parameter 


Value 


k 2 


1.3 


h 


1.9 




3.2 


v 2 


1.58 


v 3 


5 




2.5 


v s 


0.76 


h 


0.5 


k s 


0.38 


Vd 


0.95 


kd 


0.2 


n 


4 


K x 


2 


K 2 


2 


K 3 


2 


K A 


2 


Kj 


1 




0.65 



With these parameters, there are limit cycle oscillations. We leave all fixed 
except v s , and show that there are no oscillations if v s = 0.4, but oscillations 
exist if v s = 0.5 and there are delays in the negative regulatory loop, either in 
transcription or in translation (or in both). 

We choose to view the system as the feedback interconnection of two sub- 
systems, see Figure |3 
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Figure 2: Systems in feedback 



mRNA System 

The first (M) subsystem is described by the scalar differential equation 

M = v s K^/{K\ l + u 1 l)-v m M/{k m + M) 
with input u\ and output y\ — k s M. 

PER System 

The second (P) subsystem is four-dimensional: 
Po = u 2 -V 1 P /(K 1 + P a ) + V 2 P 1 /(K 2 + P 1 ) 

Pi = V 1 P Q /{K 1 + P Q )-V 2 P 1 /{K 2 + P 1 )-V 3 P 1 /(K 3 + P 1 ) + V i P 2 /(K 4 + P 2 ) 
P 2 = V 3 P 1 /(K 3 + P 1 )-V 4 P 2 /(K 4 + P 2 )-k 1 P 2 + k 2 P N -v d P 2 /(k d + P 2 ) 
Pn = hP 2 -k 2 P N 

with input u 2 and output y 2 — Pn- 
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Assume from now on that: 

v s < 0.54 (1) 

(the remaining parameters will be constrained below, in such a manner that 
those in the previously given table will satisfy all the constraints). 

As state-space for the first system, we will pick a compact interval X\ — 
[0,M], where 

< M < ^ (2) 
and we assume that v s < v m . Note that the first inequality implies that 

and therefore 

v s K?/(K?+u?) - v m M/(k m + M) < 

for all iti > 0, so that indeed X\ is forward-invariant for the dynamics. With 
the parameters shown in the table given earlier (except for v s , which is picked 
as in Q), 

M = 2.45 

satisfies all the constraints. As input space for the mRNA system, we pick 
U\ = M.>o, and as output space Y\ = [0,z; s ). Note that yi = k s M < k s M < v s . 
by (J2J), so the output belongs to Y\. 

For the second system, the state space is K> , the input space is XJi = Y\, 
and the output space is Yi = U\. 

When looking at the first system, we view U\ as ordered by the cone K.<o, 
but U<2, Yi, Y2 are all ordered in the usual manner (cone R>o). 



3 Monotonicity and Characteristics 

The first system is monotone, and has a well-defined characteristic, in the sense 
of PJ. Monotonicity is clear (one-dimensional system), and the existence of 
characteristics is immediate from the fact that M > for M < ki(ui) and 
M < for M > ki(u\), where, for each constant input u\, 

, » \ v s Kj" k m 

l[Ul) ~ v m K? + v mUl n -v,K? 

(which is an element of X\ ) . 

Note that all solutions of the differential equations which describe the M- 
system, even those that do not start in X±, enter X\ in finite time (because 
M(t) < whenever M(t) > M, for any input Ui(-)). The restriction to the 
state space Xi (instead of using all of K>o) is done for convenience, so that one 
can view the output of the M system as in input to the P-subsystem. (Desirable 
properties of the P-subsystem depend on the restriction imposed on L^-) Given 
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any trajectory, its asymptotic behavior is independent on the behavior in an 
initial finite time interval, so this does not change the conclusions to be drawn. 
(Note that solutions are defined for all times -no finite explosion times- because 
the right-hand sides of the equations have linear growth.) 

Monotonicity of the second system is also clear, from the fact that §^ > 
for all i 7^ j\ in fact, this is a strongly monotone tridiagonal system ([HIE])- 
We show that (for the parameters in the table, as well as for a larger set of 
parameters) the system has, for each constant input u, a unique equilibrium, 
and trajectories are all bounded; it follows then from [H] [7] that the unique 
equilibrium is globally asymptotically stable, which means that characteristics 
are well-defined. 

Proposition 3.1 Suppose that the following conditions hold: 

• vd + V 2 < Vi 

• Vi + Va < v 2 + v 3 

• < c < Vd 

• V 4 + v d < V3 

and that all constants are positive and the input U2(t) = c. Then the P-system 
has a unique globally asymptotically stable equilibrium. 

This will be a corollary of the following more general result. 

Theorem 1 Consider a system of the following form: 

x = c-a {x ) +f3 (xi) 
±i = a Q (x Q ) - Po(xi) - ai(xi) + f3 1 (x 2 ) 
x 2 = ax{xx) - Pi{x 2 ) -a 2 (x 2 ) -72(^2) +73(2:3) 
±3 = 72(^2) - 73(2:3) 
evolving on M> , where c > is a constant, and the functions 

«i,/3i,7i : [0,oo) -> [0,oo) 

are all differentiate, with derivatives everywhere positive, and so that oti and (3i 
are bounded, for each i, and 71,72 are unbounded. Furthermore, suppose that 
the following conditions hold: 

a 2 (oo) + /3 (oo) < a (oo) (4) 

a (oo) + /3i(oo) < ai(oo) + /3 (oo) (5) 

0:2(00) + /3i(oo) < ai(oo) (6) 

c<a 2 (oo). (7) 

Then, there is a (unique) globally asymptotically stable equilibrium for the sys- 
tem. 
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Note that (JU and J7J) imply also: 



c + /? (oo) < a (oo) . (8) 

Proof. We start by noticing that solutions are defined for all t > 0. Consider 
any maximal solution x(t) = (xo(i), x%(t), X2(t), Xs(t)). From 

— (x + xi + x 2 + x 3 ) = c-a 2 (x 2 ) (9) 

we conclude there is an estimate Xi (t) < J2 i Xi (i) < J2i x i (0) + ^ c an d hence 
there are no finite escape times. Moreover, we claim that x(-) is bounded. 

Since the system is a strongly monotone tridiagonal system, we know that 
x 3 (t) is eventually monotone. That is, for some T > 0, either 

x 3 (t) > Vi > T (10) 



or 

x 3 (t) < V* > T. (11) 

Hence, x 3 (t) admits a limit, either finite or infinite. Assume first that x 3 (t) — > 
00. Then, case (|llfl cannot hold, so l|l(J|l holds. Looking at the differential 
equation for x 3 , we know that 72 (#2 (t)) — 73( cc 3(^)) ^ f° r au t > T, which 
means that 

x 2 (t) > J2\l3( x 3 (t))) -» 00 . 

Looking again at ©, and using that c — a 2 (oo) < (property 0), we conclude 
that ^ (xo + xi + X2 + X3) (t) < for all t sufficiently large. Thus x$ + x\ + 
X2 + X3 is bounded (and nonnegative), and this implies that X2 is bounded, a 
contradiction. So x 3 is bounded. 

Next we examine the equation for X2- The two positive terms are bounded: 
the one involving a\ because a\ is a bounded function, and the one involving 
73 because X3 is bounded. Thus 

x 2 < v(t) - a 2 (x 2 ) , 

where < v(t) < k for some constant k. Thus x 2 (t) < whenever X2(t) > 
7 2 ' 1 (fc), and this proves that X2 is bounded, as claimed. 

Now we show that Xo and x\ are bounded as well. For xo, it is enough to 
notice that io < c — a (xo) + /?o(oo), so that 

x (t) > a ( 7 1 (c + /3 (oo)) =$> x o (t)<0 

so (jSJ shows that xo is bounded. Similarly, for x\ we have that xi < ao(oo) — 
Po( x i) — a i(xi) + /3i(oo) so © provides boundedness. 

Once that boundedness has been established, if we also show that there is 
a unique equilibrium then the theory of strongly monotone tridiagonal systems 
(El ED wm ensure global asymtotic stability of the equilibrium. So we show 
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that equilibria exist and are unique. It is convenient to change variables are 
write 

y Q := x + x% + x 2 + x 3 , y\ := x x + x 2 + x 3 , y 2 := x 2 + x 3 , y 3 := x 3 . 

In terms of these variables, we may set yi — 0, i = 0, 2, 1, 3, so that the equilibria 
are precisely the solutions of: 

a 2 {x 2 ) = c 

ai(xi) = a 2 (x 2 ) + (3i(x 2 ) 

a (x ) = a 2 (x 2 ) + @o(xi) 

73(^3) = 72(^2) • 

This shows uniqueness (all the functions are strictly increasing), and existence 
follows from, respectively, (7J, ©, 10} , and the fact that 73 is unbounded. I 

4 Closing the Loop 

Now we are ready to apply the main theorem in 0. In order to do this, we 
need to plot the characteristics. See Figure for the "spiderweb diagram" (the 
dotted and dashed curves are the characteristics) that shows convergence of 
the discrete iteration described in |2] when we pick the parameter v s = 0.4. 
The theorem implies that no oscillations can happen in that case, even under 
arbitrary delays in the feedback from Pn to M. 
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Figure 3: Stability of spiderweb (v s = 0.4) 

On the other hand, for a larger value, such as v s = 0.5, the discrete iteration 
conditions are violated; see Figure 0] for the "spiderweb diagram" that shows 
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Figure 4: Instability of spiderweb (v a = 0.5) 



divergence of the discrete iteration. Thus, and one may expect periodic orbits in 
this case. Indeed, simulations show that, for large enough delays, such periodic 
orbits arise, see Figure 
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Figure 5: Oscillations seen in simulations (v s = 0.5, delay of 100, initial condi- 
tions all at 0.2), using MATLAB's dde23 package 
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